
;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

begin
;
; The WRF ARW input file.  

  if (.not. isvar("input_file")) then
      start_date = "20120822" + "08"
      input_file = systemfunc("ls /public/users/wrf/wrfout/" + start_date + "/wrfout_d01_*")
  end if

  if (.not. isvar("output_prefix")) then
      output_path = "/public/users/wrf/wrfpost/test/sswdn/"
      system("mkdir -p " + output_path)
      output_prefix = output_path + "sswdn_d01_" + start_date + "_"
  end if

  a = addfile(input_file + ".nc", "r")

  wks_type = "png"
  wks_type@wkWidth = 800
  wks_type@wkHeight = 800

; Set some basic resources
  res = True
  res@MainTitle                   = "WRF SURFACE"
  res@Footer                      = False
  res@NoHeaderFooter              = True

  pltres = True
  pltres@NoHeaderFooter       = True
  mpres = True
  mpres@NoHeaderFooter       = True
  mpres@mpDataBaseVersion       = "MediumRes"
  mpres@mpDataSetName           = "Earth..4"
  mpres@mpGeophysicalLineColor      = "Black"
  mpres@mpNationalLineColor         = "Black"
  mpres@mpUSStateLineColor          = "Black"
  mpres@mpGridLineColor             = "Black"
  mpres@mpLimbLineColor             = "Black"
  mpres@mpPerimLineColor            = "Black"
  mpres@mpGeophysicalLineThicknessF = 2.0
  mpres@mpGridLineThicknessF        = 1.0
  mpres@mpLimbLineThicknessF        = 1.0
  mpres@mpNationalLineThicknessF    = 2.0
  mpres@mpUSStateLineThicknessF     = 1.0
  mpres@mpOutlineBoundarySets       = (/"National","Continental", "Geophysical"/)   ; turn on country boundaries
  mpres@mpOutlineSpecifiers = "China:states"
  mpres@mpGeophysicalLineColor      = "Black"       ; color of cont. outlines
  mpres@mpGeophysicalLineThicknessF = 2          ; thickness of outlines
  mpres@mpNationalLineThicknessF = 2          ; thickness of outlines




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Reading times
    char_times = a->Times
    hours = WRF_Times2Udunits_c(char_times, 0)
    ntimes = dimsizes(hours)         ; number of times in the file
    wrf_dt = (hours(ntimes-1) - hours(0))/(ntimes-1)
    ; convert to Beijing Time
    hours = hours + 8
    times_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")
    timesTideStr_cst = ut_string(hours, "%Y%N%D%H%M%S")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    contour_dt = 1 ; time interval of contour
    dt = tointeger(contour_dt / wrf_dt)
    do it = 0, ntimes-1, dt             ; TIME LOOP
    
    output_file = output_prefix + timesTideStr_cst(it)
    wks = gsn_open_wks(wks_type, output_file)

    print("Working on time: " + times_cst(it) )
    res@TimeLabel = times_cst(it)   ; Set Valid time to use on plots


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        

    sswdn = wrf_user_getvar(a,"SSWDN",it)    ; sswdn
      wrf_smooth_2d( sswdn, 3 )            ; smooth sswdn

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;



    ; Plotting options for SSWDN
      opts = res         
      opts@cnFillOn = True 
      opts@gsnSpreadColorEnd = -3 ; End third from the last color in color map
      opts@NoHeaderFooter       = True
      opts@cnInfoLabelOn        = False
      opts@cnLineColor = "Blue"
      opts@cnHighLabelsOn = True
      opts@cnLowLabelsOn = True
      opts@ContourParameters = (/ 0., 1100., 100. /)
      opts@cnLineLabelBackgroundColor = -1
      opts@gsnContourLineThicknessesScale = 1.5
      contour_sswdn = wrf_contour(a,wks,sswdn,opts)
      delete(opts)

    ; Contour Wind Speed
;      opts = res
;      opts@cnFillOn = True 
;      opts@ContourParameters = (/ 0., 20., 2./) 
;      opts@cnLineLabelBackgroundColor = -1
;      opts@gsnSpreadColorEnd = -3  ; End third from the last color in color map
;      contour_windspeed = wrf_contour(a,wks,wind_speed,opts)
;      delete(opts)



    ; MAKE PLOTS                                       
      plot = wrf_map_overlays(a,wks,(/contour_sswdn/),pltres,mpres)

    delete(wks)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  end do        ; END OF TIME LOOP

end
